bnlearn | Python library for learning the graphical structure | Machine Learning library

You have several errors in your CPT definitions. Primarily, you need to make sure that:

• the number of probabilities supplied are equal to the product of the number of states in the child and parent nodes,
• that the number of dimensions of the matrix/array is equal to the number of parent nodes plus one, for the child node,
• the child node should be given in the first dimension when the node dimension is greater than one.
• the names given in the dimnames arguments (e.g. the names in dimnames=list(ThisName = ...)) should match the names that were defined in the DAG, in your case with modelstring and in my answer with model2network. (So my earlier suggestion of using dimnames=list(cptNBLW = ...) should be dimnames=list(nblw = ...) to match how node nblw was declared in the model string)

You also did not add node f into your cpt list.

Below is your code with comments where things have been changed. (I have commented out the offending lines and added ones straight after)

It will be easier to let existing packages do this for you; like bnlearn. You can use custom.fit to specify the DAG and the CPTs and then use rbn to draw samples from it.

An example

How about using cpdist to draw samples from the posterior given the evidence. You can then estimate the updated parameters using bn.fit using the cpdist samples. Then plot as before.

An example:

You can add domain knowledge or constraints to structure learning in a couple of ways.

• If you want to specify the network structure and parameters using domain knowledge, you can build the network manually using custom.fit.

• If you want to estimate the structure of the BN from data then you can impose constraints on edge direction & edge presence using the whitelist and blacklist parameters in the structure learning algorithms.

• A prior can be placed on the edges in structure learning (e.g. prior="cs", where "If prior is cs, beta is a data frame with columns from, to and prob specifying the prior probability for a set of arcs. A uniform probability distribution is assumed for the remaining arcs."). There are other priors that can be used.

Regarding (1), modify takes a list or an atomic vector. bn_fit is of class bn.fit, bn.fit.dnet, however, under the hood it is a list too, since calling typeof() yields list. My guess is that there is no S3 generic method for subsetting these classes so R figures out that it is essentially a list and accordingly strips the class arguments. So subsetting bn_fit turns it into class list, and therefore you can use modify on it. Subsetting can even be done with empty brackets [], it will just return the object, but this time as class list. An alternative that I use below is to "manually" set the class attribute to NULL via attr(bnfit, "class") <- NULL.

Regarding (2), I wrote a tidyverse based function which can be used to alter the prob table of each node into a bayesm::rdirichlet distribution (see code below). The user still needs to provide part of the alpha argument (the length argument is given by the length of each prob table). Under the hood the function is relying on purrr::modify. It takes care of the classes by stripping them first and adding them back once the modification is done. My approach is to turn the prob tables into data.frames then modify the Freq column and adjust it for existing other variables (groups) and then translate that data.frame back into a table using xtabs and formulation notation via reformulate.

I'm not so deep into bayesian networks, so I do not know to what extent this function can be generalized, or whether it only works on the dataset you provided. Further, please test if the modified object is accepted by functions expecting class bn.fit, bn.fit.dnet.

I tried to comment each step of my code, but please ask if something is not clear.

(3) Regarding your question why NROW(.x) does not work in your code and you have to use NROW(node$prob) instead: This has to do with the way modify loops over the prob tables. A good way to check over what elements modify is looping over is to use purrr::pluck.

cpquery is quite difficult to work with programmatically. If you look at the examples in the help page you can see the author uses eval(parse(...)) to build the queries. I have added two approaches below, one using the methods from the help page and one using cpdist to draw samples and reweighting to get the probabilities.

Your example

With the mmhc algorithm that is used as default, you can use the layer.struct parameter to specify which pairs of layers are allowed to have edges between them. layer.struct takes a binary matrix where cell i,j is 1 if there can be edges going from variables in layer i to variables in layer j, and 0 otherwise.

The best way to use this is to combine it with the manually-specified layering of your first solution.